Spectral Graphs & Poly2Graph — A Hands‑On Tutorial

Author

Yiyang Liu

Published

October 17, 2025

0.1 Learning objectives

By the end of this tutorial, you will be able to:

  • Explain what a spectral graph is (a graph traced by the energy spectrum of a 1‑D crystal on the complex plane under open boundaries).
  • Understand the key mathematical objects that appear in modern treatments: the characteristic polynomial \(P(z,E)\), the spectral potential \(\Phi(E)\), and the density of states \(\rho(E)\).
  • Describe prior pain points (manual inspection, numerical bottlenecks, fragile extraction) and the need for an automated tool.
  • Use Poly2Graph to go from a symbolic Hamiltonian/characteristic polynomial to its spectral graph, and analyze/visualize the result.

0.1.1 Session 1

0.2 What is a spectral graph (and why is it interesting)?

Basic idea :
Consider a simple 1‑D crystal (a line of repeating “cells”). Its energy levels \(E\) depend on a wavelike parameter \(k\) (momentum).
Under periodic boundaries (“a loop”), energies trace closed loops as \(k\) runs from \(0\) to \(2\pi\).
Under open boundaries (“a finite chain”), something shocking can happen in non‑Hermitian systems (with loss/gain or non‑reciprocal couplings): the eigen‑energies collapse into arcs and loops drawn directly on the complex \(E\)‑plane. Connect those arcs/loops and you get a spectral graph.

A few key notions :

  • The physics lives in the characteristic polynomial of the Bloch Hamiltonian \(H(z)\):
    \(P(z,E)=\det[H(z)-E I] = \sum_{n=-p}^q a_n(E)\,z^n, \quad z=e^{ik}\)
  • Sort the \(z\)‑roots \(\beta_i(E)\) of \(P(z,E)=0\) by magnitude. On the spectral graph, the “middle two” obey \(|\beta_p(E)|=|\beta_{p+1}(E)|\). This condition defines the generalized Brillouin zone (GBZ) and pins down the spectral loci under OBC.
  • Define a spectral potential (a Ronkin function) from the roots (one convenient form):
    \(\Phi(E) = - \log|a_q(E)| - \sum_{i=p+1}^{p+q}\log|\beta_i(E)|\) Its Laplacian gives the density of states (DOS):
    \(\rho(E)= -\tfrac{1}{2\pi}\nabla^2\Phi(E)\)

Intuition: treat each eigenvalue \(\epsilon_n\) as a tiny electric charge \(\rho(E) = \lim_{N\to\infty}\sum_n \frac{1}{N} \delta(E-\epsilon_n)\); while \(\Phi(E)\) is an electrostatic‑like potential whose “ridges/troughs” trace the spectral graph; \(\rho(E)\) is nonzero only on the graph. Please see the Section 1.0.1 below for a concrete illustration.

Why interesting?

Non‑Hermitian spectra can realize many more graph connectivities than standard topological labels (like \(\mathbb{Z},\mathbb{Z}_2\)): stars, flowers, braids, etc., with novel transitions that change the number of endpoints/junctions/loops (no Hermitian analog). This is naturally expressed by writing dispersions as bivariate Laurent polynomials \(P(E,z)\) and studying how their imaginary‑momentum (“inverse skin depth”) surfaces intersect across the \(E\)‑plane.
Moreover, spectral graph is a general algebra-to-graph bridge. Because the equation \(\Phi(E) = - \log|a_q(E)| - \sum_{i=p+1}^{p+q}\log|\beta_i(E)|\) uses only \(P(E,z)\), the construction extends to arbitrary Laurent polynomials, and via standard reductions to vectors and matrices. You get a universal topological fingerprint for many algebraic objects—making graph learning a front-door to problems in linear algebra, signal processing, photonics, circuits, and more.
Here, we visualize a single-band example by an interactive 3D surface plot of \(\Phi(E)\) along with the extracted skeleton Figure 1.

0.2.1 Session 2

0.3 What was hard before (and why we need Poly2Graph)?

Prior to Poly2Graph, extracting spectral graphs was manual and fragile: - You had to scan parameters, plot large finite‑chain spectra, and visually infer the continuum graph. - Computing the GBZ roots \(\beta_i(E)\) for every sampled energy \(E\) is a bottleneck; naive solvers are slow and memory‑hungry. - Even with a DOS image, turning that into a clean graph with correct endpoints, junctions, and loops requires careful morphology (thresholding, skeletonization, node/edge parsing).

Poly2Graph lifts these barriers with an end‑to‑end, high‑performance pipeline: 1. Accepts either \(H(z)\) (as a sympy.Matrix) or \(P(z,E)\) (string or sympy.Poly).
2. Computes the spectral potential \(\Phi(E)\) from \(z\)‑roots (fast eigen‑solvers on Frobenius companion matrices; optional GPU/TPU via TensorFlow/PyTorch).
3. Takes the Laplacian to get DOS, then does adaptive, two‑stage refinement (coarse mask → high‑res on just ~1–5% of the region).
4. Applies image morphology to extract a one‑pixel skeleton, then returns a NetworkX MultiGraph with full edge geometry (the sampled curve).
This automation has enabled HSG‑12M (11.6M static + 5.1M temporal spectral multigraphs), demonstrating scale and diversity previously out of reach(Learn more from @arXiv:2506.08618 if interested).

Now let’s try to use Poly2Graph ***

# If running on devices without essential libraries, uncomment the next line:
# !pip install -U poly2graph sympy networkx matplotlib

import sympy as sp
import numpy as np
np.set_printoptions(threshold=20)
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import networkx as nx

# Poly2Graph (must be installed in your environment)
# If this import fails, run the pip command above.
import poly2graph as p2g
# For better plotly rendering in Jupyter environments
import plotly.io as pio
pio.renderers.default = "plotly_mimetype+notebook_connected"

Please ensure your version of Python >= 3.11 , you can check the installation:

print("Poly2Graph version:", p2g.__version__)
Poly2Graph version: 0.2.0
# Always remember to define common symbols first
k = sp.symbols('k', real=True)
z, E = sp.symbols('z E', complex=True)

0.3.1 1. Quick start

0.4 one‑band example

We’ll start from a simple one‑band characteristic polynomial and obtain its spectral graph.

Characteristic polynomial $ P(z,E) = z^4 - z^2 - z^{-1} - E $ (This corresponds to a scalar Bloch Hamiltonian \(h(z)=z^4-z^2-z^{-1}\), with \(z=e^{ik}\).)

The valid input formats to initialize a p2g.SpectralGraph object are: 1. Characteristic polynomial in terms of z and E: - as a string of the Poly in terms of z and E - as a sympy’s Poly (sympy.polys.polytools.Poly) with {z, 1/z, E} as generators 2. Bloch Hamiltonian in terms of k or z - as a sympy Matrix in terms of k - as a sympy Matrix in terms of z

All the following characteristics are valid and will initialize to the same characteristic polynomial and therefore produce the same spectral graph

char_poly_str = '-z**-1 - E - z**2 + z**4'

char_poly_Poly = sp.Poly(
    -z**-1 - E - z**2 + z**4,
    z, 1/z, E # generators are z, 1/z, E
)

phase_k = sp.exp(sp.I*k)
char_hamil_k = sp.Matrix([-phase_k**-1 - phase_k**2 + phase_k**4])

char_hamil_z = sp.Matrix([-z**-1 - E - z**2 + z**4])

Here we use the string format as the initialization example.

sg = p2g.SpectralGraph(char_poly_str, k=k, z=z, E=E)
# print the number of energy bands
print('Bands:', sg.num_bands)

# print the max hopping range (p,q)
print('Max hop (right):', sg.poly_p, 'Max hop (left):', sg.poly_q)

# print the [Ereal_min,Ereal_max,Eimag_min,Eimag_max] window
print('Spectral window (auto-estimated):', sg.spectral_square) 
Bands: 1
Max hop (right): 1 Max hop (left): 4
Spectral window (auto-estimated): [-2.31139048  1.6398687  -1.97562959  1.97562959]

We will see Poly2Graph automatically compute set of properties next:

sg.ChP

\(\displaystyle \operatorname{Poly}{\left( z^{4} - z^{2} - \frac{1}{z} - E, z, \frac{1}{z}, E, domain=\mathbb{Z} \right)}\)

Bloch Hamiltonian:

  • For one-band model(The exponent of \(E\) is 1), it is a unique, rank-0 matrix (scalar)
sg.h_k

\(\displaystyle \left[\begin{matrix}e^{4 i k} - e^{2 i k} - e^{- i k}\end{matrix}\right]\)

or in \(z\) format(recall that \(z = e^{ik}\))

sg.h_z

\(\displaystyle \left[\begin{matrix}- \frac{- z^{5} + z^{3} + 1}{z}\end{matrix}\right]\)

Frobenius Companion Matrix of \(P(z, E)\):

  • The Frobenius companion matrix is a special matrix whose eigenvalues are exactly the roots of the characteristic polynomial \(P(z, E)\) for fixed \(E\).
  • In Poly2Graph, sg.companion_E gives this matrix (with \(E\) as a parameter, \(z\) as variable).
  • Purpose:
    • Efficiently computes all \(z\)-roots for any \(E\) using linear algebra (eigenvalues), instead of slow root-finding.
    • These \(z\)-roots are essential for calculating the spectral graph, GBZ, spectral potential, and density of states.
sg.companion_E

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0 & 1\\1 & 0 & 0 & 0 & E\\0 & 1 & 0 & 0 & 0\\0 & 0 & 1 & 0 & 1\\0 & 0 & 0 & 1 & 0\end{matrix}\right]\)

It is a key tool for automating and accelerating spectral graph extraction in Poly2Graph.

0.4.1 Spectral potential \(\Phi(E)\), density of states \(\rho(E)\), and skeleton

We compute \(\Phi(E)\) and \(\rho(E)\) on a grid of complex energies, then binarize \(\rho(E)\) to get a skeleton.

::: {#cell-1 band spectral images .cell execution_count=12}

phi, dos, skel = sg.spectral_images()

fig, axes = plt.subplots(1, 3, figsize=(8, 3), sharex=True, sharey=True)

# 1) Spectral potential
axes[0].imshow(phi, extent=sg.spectral_square, origin='lower', cmap='terrain')
axes[0].set(xlabel='Re(E)', ylabel='Im(E)', title='Spectral Potential')

# 2) Density of States
pmin, pmax = np.percentile(dos, (0.1, 99.9))
# Clip extreme DOS to increase visibility
norm = colors.Normalize(vmin=pmin, vmax=pmax)
axes[1].imshow(dos, extent=sg.spectral_square, cmap='viridis', norm=norm)
axes[1].set(xlabel='Re(E)', title='Density of States')

# 3) Binarized skeleton
axes[2].imshow(skel, extent=sg.spectral_square, cmap='gray')
axes[2].set(xlabel='Re(E)', title='Graph Skeleton')
plt.tight_layout()
plt.savefig("assets/1-band-example_si.png")
plt.show()

Spectral images for 1 band example: (1) Spectral potential Φ, (2) Density of States ρ, (3) Graph skeleton from ρ

:::

0.4.2 Extract the spectral graph

Next, we visualize these results by a interactive 3D surface plot of \(\Phi(E)\) along with the extracted skeleton.

import plotly.graph_objects as go

E_re = np.linspace(*sg.spectral_square[:2], phi.shape[0])
E_im = np.linspace(*sg.spectral_square[2:], phi.shape[1])
fig = go.Figure(data=[go.Surface(z=phi, x=E_re, y=E_im, 
            opacity=0.6, colorscale='Spectral_r')])
fig.update_layout(
scene=dict(
    aspectratio=dict(x=1.5, y=1, z=.8),
    xaxis_title='Re(E)',yaxis_title='Im(E)',zaxis_title='Phi(E)',
),
title='3D Surface Plot of Spectral Potential',
width=700, height=600
)
fig.show()
Figure 1: 3D surface plot

sg.spectral_graph() returns a NetworkX MultiGraph with node positions and polyline geometry along each edge.

::: {#cell-1 band spectral graph .cell execution_count=14}

# Draw the multi-band spectral graph
G = sg.spectral_graph()
pos2 = nx.get_node_attributes(G, 'pos')
plt.figure()
nx.draw_networkx_nodes(G, pos2, node_size=20)
nx.draw_networkx_edges(G, pos2)
plt.axis('equal'); plt.title('Spectral Graph (1-band)')
plt.show()

print(f'#nodes: {G.number_of_nodes()}  #edges: {G.number_of_edges()}')

# Show attributes of a node and an edge
any_node = list(G.nodes(data=True))[0]
print('A node has attributes:', list(any_node[1].keys()))
any_edge = list(G.edges(keys=True, data=True))[0]
print('An edge has attributes:', list(any_edge[3].keys()))

Spectral Graph (1-band)
#nodes: 7  #edges: 6
A node has attributes: ['pos', 'dos', 'potential']
An edge has attributes: ['weight', 'pts', 'avg_dos', 'avg_potential']

:::

As you can see above, the spectral graph is a networkx.MultiGraph object.

  • Node Attributes
    1. pos : (2,)-numpy array
      • the position of the node \((\text{Re}(E), \text{Im}(E))\)
    2. dos : float
      • the density of states at the node
    3. potential : float
      • the spectral potential at the node
  • Edge Attributes
    1. weight : float
      • the weight of the edge, which is the length of the edge in the complex energy plane
    2. pts : (w, 2)-numpy array
      • the positions of the points constituting the edge, where w is the number of points along the edge, i.e., the length of the edge, equals weight
    3. avg_dos : float
      • the average density of states along the edge
    4. avg_potential : float
      • the average spectral potential along the edge
node_attr = dict(G.nodes(data=True))
edge_attr = list(G.edges(data=True))
print('The attributes of the first node\n', node_attr[0], '\n')
print('The attributes of the first edge\n', edge_attr[0][-1], '\n')
The attributes of the first node
 {'pos': array([0.88936098, 1.90810319]), 'dos': np.float32(0.25110298), 'potential': np.float32(-0.4699422)} 

The attributes of the first edge
 {'weight': np.float64(2.0950480742544872), 'pts': array([[0.88550233, 1.89652723],
       [0.88550233, 1.89266858],
       [0.88550233, 1.88880993],
       ...,
       [0.81990525, 0.01350528],
       [0.8237639 , 0.00964663],
       [0.8237639 , 0.00192933]]), 'avg_dos': np.float32(0.079228334), 'avg_potential': np.float32(-0.2639912)}

0.4.3 2. Multi‑band example

0.5 Multi‑band example

Now consider a four‑band polynomial:

\[ P(z,E) = z^2 + z^{-2} + E z - E^4 \]

A possible Bloch Hamiltonian realization is: \[\textbf{h}(z)=\begin{pmatrix} 0 & 0 & 0 & z^2 + 1/z^2 \\ 1 & 0 & 0 & z \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ \end{pmatrix}\] Poly2Graph can work directly from the string/sympy.Poly as shown below.

# We start with a polynomial string
sg_multi = p2g.SpectralGraph("z**2 + 1/z**2 + E*z - E**4", k, z, E)
# print the number of energy bands
print('Bands:', sg_multi.num_bands)

# print the max hopping range (p,q)
print('Max hop (right):', sg_multi.poly_p, 'Max hop (left):', sg_multi.poly_q)

# print the [Ereal_min,Ereal_max,Eimag_min,Eimag_max] window
print('Spectral window (auto-estimated):', sg_multi.spectral_square)
Bands: 4
Max hop (right): 2 Max hop (left): 2
Spectral window (auto-estimated): [-1.40559546  1.40559546 -1.40559546  1.40559546]
sg_multi.ChP

\(\displaystyle \operatorname{Poly}{\left( z^{2} + zE + \frac{1}{z^{2}} - E^{4}, z, \frac{1}{z}, E, domain=\mathbb{Z} \right)}\)

sg_multi.h_k

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & e^{2 i k} + e^{- 2 i k}\\1 & 0 & 0 & e^{i k}\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\end{matrix}\right]\)

sg_multi.h_z

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & z^{2} + \frac{1}{z^{2}}\\1 & 0 & 0 & z\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\end{matrix}\right]\)

sg_multi.companion_E

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & -1\\1 & 0 & 0 & 0\\0 & 1 & 0 & E^{4}\\0 & 0 & 1 & - E\end{matrix}\right]\)

0.5.1 Spectral potential \(\Phi(E)\), density of states \(\rho(E)\), and skeleton

Like the one-band example, we also compute \(\Phi(E)\) and \(\rho(E)\) on a grid of complex energies, then binarize \(\rho(E)\) to get a skeleton.

::: {#cell-Multi-band spectral images .cell execution_count=22}

phi2, dos2, skel2 = sg_multi.spectral_images()

fig, axes = plt.subplots(1, 3, figsize=(8, 3), sharex=True, sharey=True)

# 1) Spectral potential
axes[0].imshow(phi2, extent=sg.spectral_square, origin='lower', cmap='terrain')
axes[0].set(xlabel='Re(E)', ylabel='Im(E)', title='Spectral Potential')

# 2) Density of States
pmin, pmax = np.percentile(dos, (0.1, 99.9))
# Clip extreme DOS to increase visibility
norm = colors.Normalize(vmin=pmin, vmax=pmax)
axes[1].imshow(dos2, extent=sg.spectral_square, cmap='viridis', norm=norm)
axes[1].set(xlabel='Re(E)', title='Density of States')

# 3) Binarized skeleton
axes[2].imshow(skel2, extent=sg.spectral_square, cmap='gray')
axes[2].set(xlabel='Re(E)', title='Graph Skeleton')
plt.tight_layout()
plt.savefig("assets/multi-band-example_si.png")
plt.show()

Spectral images for multi-band example: (1) Spectral potential Φ, (2) Density of States ρ, (3) Graph skeleton from ρ

:::

0.5.2 Extract the spectral graph

::: {#cell-Multi-band spectral graph .cell execution_count=23}

# Draw the multi-band spectral graph
G2 = sg_multi.spectral_graph()
pos2 = nx.get_node_attributes(G2, 'pos')
plt.figure()
nx.draw_networkx_nodes(G2, pos2, node_size=20)
nx.draw_networkx_edges(G2, pos2)
plt.axis('equal'); plt.title('Spectral Graph (multi‑band)')
plt.show()

print(f'#nodes: {G2.number_of_nodes()}  #edges: {G2.number_of_edges()}')

# Show attributes of a node and an edge
any_node = list(G2.nodes(data=True))[0]
print('A node has attributes:', list(any_node[1].keys()))
any_edge = list(G2.edges(keys=True, data=True))[0]
print('An edge has attributes:', list(any_edge[3].keys()))

Spectral Graph (multi‑band)
#nodes: 21  #edges: 18
A node has attributes: ['pos', 'dos', 'potential']
An edge has attributes: ['weight', 'pts', 'avg_dos', 'avg_potential']

:::

0.5.3 3. Characteristic Polynomial Class

0.6 p2g.CharPolyClass

For scanning families of polynomials (e.g., varying two complex coefficients across a grid), use p2g.CharPolyClass. This class parallelizes evaluation and is designed to generate many graphs efficiently (e.g., to reproduce temporal sequences as in T‑HSG‑5M).

To be specific, p2g.CharPolyClass is: - A class of parametrized non-Hermitian lattices - Generate spectral properties, including spectral graph in parallel - Optimized for computational efficiency and numerical stability

Let us add two parameters {a,b} to the aforementioned multi-band example and construct a p2g.CharPolyClass object.

a, b = sp.symbols('a b', real=True)

cp = p2g.CharPolyClass(
    "z**2 + a/z**2 + b*E*z - E**4", 
    k=k, z=z, E=E,
    params={a, b}, # pass parameters as a set
)
Derived Bloch Hamiltonian `h_z` with 4 bands.

Just like the examples above, we can also check the automatically computed properties:

cp.ChP

\(\displaystyle \operatorname{Poly}{\left( z^{2} + b zE + a \frac{1}{z^{2}} - E^{4}, z, \frac{1}{z}, E, domain=\mathbb{Z}\left[a, b\right] \right)}\)

cp.h_k

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & a e^{- 2 i k} + e^{2 i k}\\1 & 0 & 0 & b e^{i k}\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\end{matrix}\right]\)

cp.h_z

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & \frac{a}{z^{2}} + z^{2}\\1 & 0 & 0 & b z\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\end{matrix}\right]\)

cp.companion_E

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & - a\\1 & 0 & 0 & 0\\0 & 1 & 0 & E^{4}\\0 & 0 & 1 & - E b\end{matrix}\right]\)

To get an array of spectral images or spectral graphs, we first prepare the values of the parameters {a,b}

a_array = np.linspace(-2, 2, 6)
b_array = np.linspace(-1, 1, 6)
a_grid, b_grid = np.meshgrid(a_array, b_array)
param_dict = {a: a_grid, b: b_grid}
print('a_grid shape:', a_grid.shape,
    '\nb_grid shape:', b_grid.shape)
a_grid shape: (6, 6) 
b_grid shape: (6, 6)

You can check the grid generated to ensure your parameters

a_grid, b_grid
(array([[-2. , -1.2, -0.4,  0.4,  1.2,  2. ],
        [-2. , -1.2, -0.4,  0.4,  1.2,  2. ],
        [-2. , -1.2, -0.4,  0.4,  1.2,  2. ],
        [-2. , -1.2, -0.4,  0.4,  1.2,  2. ],
        [-2. , -1.2, -0.4,  0.4,  1.2,  2. ],
        [-2. , -1.2, -0.4,  0.4,  1.2,  2. ]]),
 array([[-1. , -1. , -1. , -1. , -1. , -1. ],
        [-0.6, -0.6, -0.6, -0.6, -0.6, -0.6],
        [-0.2, -0.2, -0.2, -0.2, -0.2, -0.2],
        [ 0.2,  0.2,  0.2,  0.2,  0.2,  0.2],
        [ 0.6,  0.6,  0.6,  0.6,  0.6,  0.6],
        [ 1. ,  1. ,  1. ,  1. ,  1. ,  1. ]]))

Note that the value array of the parameters should have the same shape, which is also the shape of the output array of spectral images

phi_arr, dos_arr, binaried_dos_arr, spectral_square = \
    cp.spectral_images(param_dict=param_dict, device='/cpu:0')
print('phi_arr shape:', phi_arr.shape,
    '\ndos_arr shape:', dos_arr.shape,
    '\nbinaried_dos_arr shape:', binaried_dos_arr.shape)
phi_arr shape: (6, 6, 1024, 1024) 
dos_arr shape: (6, 6, 1024, 1024) 
binaried_dos_arr shape: (6, 6, 1024, 1024)
phi_arr[:, :, 0, 0]
array([[-2.79088283, -2.56152201, -2.07190657, -2.33312154, -2.826159  ,
        -3.05699062],
       [-2.61931658, -2.3819747 , -1.87466085, -2.14679098, -2.65622234,
        -2.89444447],
       [-2.42141795, -2.17299414, -1.64018846, -1.92253423, -2.45588088,
        -2.70452452],
       [-2.42141819, -2.17299461, -1.64018691, -1.92253542, -2.45588088,
        -2.70452452],
       [-2.61931658, -2.3819747 , -1.87466145, -2.14679122, -2.65622234,
        -2.8944447 ],
       [-2.79088283, -2.56152201, -2.07190657, -2.33312154, -2.82615995,
        -3.05699062]])

phi_arr, dos_arr, and binaried_dos_arr are arrays computed by cp.spectral_images(), representing different spectral quantities:

  • phi_arr: The array of spectral potential values, i.e., \(\Phi(E)\), for each parameter point((a,b), here \(6*6\)) and energy grid(each small sampling grid on the complex \(E\) plain, here \(1024*1024\)). It describes the “potential landscape” over the complex energy plane, with its ridges tracing the spectral graph.
  • dos_arr: The array of density of states values, i.e., \(\rho(E)\), for each parameter point and energy grid. It is significant only on the spectral graph (where the spectral condition is met), and nearly zero elsewhere.
  • binaried_dos_arr: The binarized density of states array, obtained by thresholding \(\rho(E)\), which extracts the skeleton (the graph structure) for further analysis.

These arrays typically have shapes like (number of parameter points, ..., H, W), where \(H\) and \(W\) are the height and width of the energy grid. Each parameter point corresponds to one spectral potential/density/skeleton image.

::: {#cell-ChP spectral images .cell execution_count=33}

from mpl_toolkits.axes_grid1 import ImageGrid

fig = plt.figure(figsize=(13, 13))
grid = ImageGrid(fig, 111, nrows_ncols=(6, 6), axes_pad=0, 
                 label_mode='L', share_all=True)

for ax, (i, j) in zip(grid, [(i, j) for i in range(6) for j in range(6)]):
    ax.imshow(phi_arr[i, j], extent=spectral_square[i, j], cmap='Spectral_r')
    ax.set(xlabel='Re(E)', ylabel='Im(E)')
    ax.text(
        0.03, 0.97, f'a = {a_array[i]:.2f}, b = {b_array[j]:.2f}',
        ha='left', va='top', transform=ax.transAxes,
        fontsize=10, color='tab:red',
        bbox=dict(alpha=0.8, facecolor='white')
    )

plt.tight_layout()
plt.savefig("assets/ChP-example_si.png")
plt.show()
C:\Users\ReYangonLiu\AppData\Local\Temp\ipykernel_6764\1050636169.py:20: UserWarning:

This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.

a & b parameter grid of spectral potential Φ(E)

:::